#-----------------------------------#
# Replication of:
# Seeking Opportunity in the Knowledge Economy: Moving Places, Moving Politics?
# Valentina Consiglio, Thomas Kurer

# File: 8_models_main_LLMR

# Note: This file re-runs the main analysis from script 5_models_main_NUTS3
# at the level of labor market regions (LMR). 
#-----------------------------------#


# Clear workspace
rm(list=ls())

#### Load Packages ####
library(pacman)
p_load(car, arm, lmtest, vars, ggplot2, psych, stargazer, dplyr, readxl, tidyverse, 
       haven, labelled, gtsummary, survey, srvyr, BMisc, questionr, plm, texreg, broom, fixest, modelsummary, bife, data.table)

#-------------------------------------------------------------------------------# 

#### Define relevant paths ####
data_path <- "data_created/"
graph_path <- "figures/" 
tables_path <- "tables/" 

#### Import data #### 

load((paste0(data_path, "df_2009_soepoi_v4.RData")))


#-----------------------------------------------------------------------------------------------------#
#### Additional data manipulation and creation of subsamples ####
#-----------------------------------------------------------------------------------------------------#

#-----------------------------------------------------------------------------------------------------#
##### Subsets of data for different years ####
#-----------------------------------------------------------------------------------------------------#

# df_2009: start year 2009

# df: standard, start year 2010
df <- df_2009 %>%
  filter(syear>2009)

# df_2014: start year 2014
df_2014 <- df_2009 %>%
  filter(syear>2013)

#-----------------------------------------------------------------------------------------------------#
##### Additional variables ####
#-----------------------------------------------------------------------------------------------------#

# Store datasets in list
datasets <- list(df_2009 = df_2009, df = df, df_2014 = df_2014)


# Apply data transformation to each dataset
datasets <- lapply(datasets, function(data) {
  data %>%
    group_by(pid) %>%
    mutate(
      # Create upmove and downmove
      upmove_ub_lmr = ifelse(opmove_ub_lmr == 1, 1,
                             ifelse(opmove_ub_lmr %in% c(2,3), 0, NA)),
      
      downmove_ub_lmr = ifelse(opmove_ub_lmr == 3, 1,
                             ifelse(opmove_ub_lmr %in% c(1,2), 0, NA)),
      # Create index counting number of observations 
      count_years = syear - first(syear) + 1,
      
      # New vars with 0 for first year and NA = NA
      upmove_ub_lmr2 = ifelse(count_years== 1, 0, # take value NA if original takes NA
                              ifelse( is.na(upmove_ub_lmr), NA, upmove_ub_lmr)), # take value 0 if first year in dataset 
      downmove_ub_lmr2 = ifelse(count_years== 1, 0,
                                ifelse(is.na(downmove_ub_lmr), NA, downmove_ub_lmr)), 
      lmr_moved_ub2 = ifelse(count_years== 1, 0, 
                             ifelse(is.na(lmr_moved_ub), NA, lmr_moved_ub)),
      
      
      # Cumulative sum (count move variable, based on variable that is 0 in first year)
      count_lmr_moved_ub2 = cumsum(ifelse(is.na(lmr_moved_ub2), 0, lmr_moved_ub2)),
      
      # "Stay treated" variable (remains 1 from first year of moving)
      lmr_moved_ub2_cont = ifelse(cumsum(ifelse(is.na(lmr_moved_ub2), 0, lmr_moved_ub2)) > 0, 1, 0),
      
      # Total up / down
      
      total_up_ub2 = sum(replace_na(upmove_ub_lmr2, 0)),
      total_down_ub2 = sum(replace_na(downmove_ub_lmr2, 0)),
      
      # Flag multiple movers
      mm_ub2 = ifelse(total_up_ub2 + total_down_ub2 > 1, 1, ifelse(total_up_ub2 + total_down_ub2 %in% c(0,1), 0, NA)),
      
      # Recode opportunity index for symmetry analysis (downmoves)
      opportunity_KE_pc01_lmr_rc = opportunity_KE_pc01_lmr * (-1)
      
    ) %>%
    ungroup()
})

# Transform datasets into original datasets
df_2009 <- datasets$df_2009
df <- datasets$df
df_2014 <- datasets$df_2014


#-----------------------------------------------------------------------------------------------------#
##### Subsamples upmove/downmove ####
#-----------------------------------------------------------------------------------------------------#

# Keep only upward moves in sample 

# -- based on unbalanced vars
datasets_upmove_ub2 <- lapply(datasets, function(data) {
  data %>%
    filter((total_up_ub2 == 1 | total_up_ub2 == 0) & mm_ub2 == 0 & total_down_ub2 == 0)
  
})

# Store datasets from list
df_2009_upmove_ub2 <- datasets_upmove_ub2$df_2009
df_upmove_ub2 <- datasets_upmove_ub2$df
df_2014_upmove_ub2 <- datasets_upmove_ub2$df_2014

rm(datasets_upmove2, datasets_upmove_ub2) # remove list with datasets


# Keep only downward moves in sample (based on balanced vars)

# -- based on unbalanced vars
datasets_downmove_ub2 <- lapply(datasets, function(data) {
  data %>%
    filter((total_down_ub2 == 1 | total_down_ub2 == 0) & mm_ub2 == 0 & total_up_ub2 == 0)
  
})

# Store datasets from list
df_2009_downmove_ub2 <- datasets_downmove_ub2$df_2009
df_downmove_ub2 <- datasets_downmove_ub2$df
df_2014_downmove_ub2 <- datasets_downmove_ub2$df_2014


rm(datasets_downmove2, datasets_downmove_ub2, datasets) # remove lists with datasets


#-----------------------------------------------------------------------------------------------------#
##### Subsamples without multiple movers ####
#-----------------------------------------------------------------------------------------------------#

## Dataset without multiple movers ##

# Unbalanced
df_womm <- df %>%
  filter(mm_ub2 == 0) # drops around 6000 indXyear observations

df_2009_womm <- df_2009 %>%
  filter(mm_ub2 == 0)

df_2014_womm <- df_2014 %>%
  filter(mm_ub2 == 0)


#-------------------------------------------------------------------------------# 
#### 1 | FE-Models with Opportunity Index at LMR-level as IV  ####
#-------------------------------------------------------------------------------# 

## Texreg input 
vec_gof_names <- c("Num. obs." = "N",
                   "Num. groups: pid" = "N individuals",
                   "Num. groups: syear" = "N years")



#-------------------------------------------------------------------------------# 
##### 1.1 | Political integration and orientation #####
#-------------------------------------------------------------------------------# 


fe_main_polint_polor2_lmr <-  list(
  "Volunteering (Yes/No)" = feols(ft_volunt_d ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) 
                                  | pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df),
  
  "Local political activism (Yes/No)" = feols(ft_localpol_d ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) 
                                              | pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df),
  
  "Voting in federal election (Yes/No)" = feols(votefedel_lead ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)
                                                | pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df_2009),
  
  "Party leaning (Yes/No)" = feols(partyleand2 ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)
                                   | pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df),
  
  "Party leaning intensity (0-5)" = feols(partyleanint ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)
                                          | pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df),
  
  "Political orientation (0 L – 10 R)" = feols(polor ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) 
                                               | pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df_2009),
  
  "Centre left party ID" = feols(centre_left2 ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) 
                                 | pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df),
  
  "Centre right party ID" = feols(centre_right2 ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) |
                                    pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df)
)


texreg(fe_main_polint_polor2_lmr)


texreg(fe_main_polint_polor2_lmr, digits = 3,
       label = "tab:fe_main_polint_polor2_lmr",
       caption = "Opportunity moves and political integration and orientation (LMR)",
       custom.coef.map = list('opportunity_KE_pc01_lmr' = 'Opportunity Index (LMR)'),
       custom.gof.rows = list("Individual fixed-effects" = c("YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES"), "Time fixed-effect" = c("YES", "YES", "YES", "YES", "YES", "YES", "YES", "YES")),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. All models include age group, education group, and household type as control variables; Standard errors are clustered at the LMR-level. Source: SOEP v.37, 2009/10-2020."),
       file = paste0(tables_path, "robustness_heterogeneity/Appendix_Table_D27.tex"))



#--------------------------------------------------------------------------------------------#

#-------------------------------------------------------------------------------# 
##### 1.2 | Party identification #####
#-------------------------------------------------------------------------------# 


# With alternative coding where all those with 0 party identification are in the reference category of the dummies

fe_main_partid2_lmr <-  list(
  "SPD " = feols(spd_co2 ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                   pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df),
  "CDU/CSU " = feols(cducsu_co2 ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                       pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df),
  "Gruene " = feols(greens_co2 ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                      pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df),
  "FDP " = feols(fdp_co2 ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                   pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df),
  "Linke" = feols(linke_co2 ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                    pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df),
  "AfD " = feols(afd_co2 ~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype) | 
                   pid + syear, panel.id = ~pid + syear, cluster = "lmr_cd", data = df_2014)
)

texreg(fe_main_partid2_lmr)


texreg(fe_main_partid2_lmr, digits = 3,
       label = "tab:fe_main_partid2_lmr",
       caption = "Opportunity moves and party identification (LMR)",
       custom.coef.map = list('opportunity_KE_pc01_lmr' = 'Opportunity Index (LMR)'),
       custom.gof.rows = list("Individual fixed-effects" = c("YES", "YES", "YES", "YES", "YES", "YES"), "Time fixed-effect" = c("YES", "YES", "YES", "YES", "YES", "YES")),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                            
       custom.note = paste("%stars. Party leanings include coalitions in the political direction of the respective party (i.e., including left coalitions for SPD, Gruene, and Linke; right coalitions for CDU/CSU, FDP, and AfD). All models include age group, education group, and household type as control variables; Standard errors are clustered at the LMR-level. Source: SOEP v.37, 2010-2020."),
       file = paste0(tables_path, "robustness_heterogeneity/Appendix_Table_D28.tex"))





#-------------------------------------------------------------------------------# 
#### 2 | Economic and Cultural mechanisms ####

#-------------------------------------------------------------------------------# 

##### 2.1 | Economic #####


# Main #
econ_DVs <- c("ln_pglabgro2", "pgisei08", "worr_wpsec")

names_econ_DVs <- c(ln_pglabgro2 = "Log gross monthly earnings", 
                    pgisei08 = "ISEI-score",
                    worr_wpsec = "Worries about workplace security"#, 
                    #income = "Equiv. HH income"
)


fe_mech_econ_lmr <- list()

# Loop over econ CVs
for (econ_DV in econ_DVs) {
  
  # Define model
  formula <- as.formula(paste(econ_DV, "~ opportunity_KE_pc01_lmr + pgexpft + pgexpft^2 + as.factor(educ_h) + as.factor(hhtype)| pid + syear"))
  
  # Run regression 
  fe_mech_econ_lmr[[names_econ_DVs[[econ_DV]]]] <- feols(formula, panel.id = ~pid + syear, cluster = "lmr_cd", data = df_2009)
}

texreg(fe_mech_econ_lmr)
texreg(fe_mech_econ_lmr, digits = 3,
       label = "tab:fe_mech_econ_lmr",
       caption = "Opportunity moves: Socio-economic (LMR)",
       custom.coef.map = list('opportunity_KE_pc01_lmr' = 'Opportunity Index (LMR)'),
       custom.gof.rows = list("Individual fixed-effects" = rep("YES", length(fe_mech_econ_lmr)), 
                              "Time fixed-effect" = rep("YES", length(fe_mech_econ_lmr))),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. All models include work experience, work experience squared, education group, and household type as control variables; Standard errors are clustered at the LMR-level. Source: SOEP v.37, 2009-2020."),
       file = paste0(tables_path, "robustness_heterogeneity/Appendix_Table_D29a.tex"))


#------------------------------------------------------------------------------#
##### 2.2 | Cultural #####
#------------------------------------------------------------------------------#

# DVs on cultural variables

cult_DVs <- c("ft_opera_d", "ft_cinema_d","ft_artmus_d") # 

names_cult_DVs <- c(ft_opera_d = "Cultural events (dummy)", 
                    ft_cinema_d = "Cinema, concerts, clubs (dummy)", 
                    ft_artmus_d = "Artistic/musical activities (dummy)"
)

# , ft_index_cult = "Cultural activities (index)"

fe_mech_culture_lmr <- list()

# Loop over cult CVs
for (cult_DV in cult_DVs) {
  
  # Define model
  formula <- as.formula(paste(cult_DV, "~ opportunity_KE_pc01_lmr + as.factor(age_grp3) + as.factor(educ_h) + as.factor(hhtype)| pid + syear"))
  
  # Run regression 
  #  model_results_cult[[cult_DV]] <- feols(formula, panel.id = ~pid + syear, cluster = "lmr_cd", data = df)
  fe_mech_culture_lmr[[names_cult_DVs[[cult_DV]]]] <- feols(formula, panel.id = ~pid + syear, cluster = "lmr_cd", data = df_2009)
}


texreg(fe_mech_culture_lmr)
texreg(fe_mech_culture_lmr, digits = 3,
       label = "tab:fe_mech_culture_lmr",
       caption = "Opportunity moves: Cultural Outcomes (LMR)",
       custom.coef.map = list('opportunity_KE_pc01_lmr' = 'Opportunity Index (LMR)'),
       custom.gof.rows = list("Individual fixed-effects" = rep("YES", length(fe_mech_culture_lmr)), 
                              "Time fixed-effect" = rep("YES", length(fe_mech_culture_lmr))),
       stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.adjrs = F,
       custom.gof.names = vec_gof_names,                                             
       custom.note = paste("%stars. All models are linear probability models and include age group, education group, and household type as control variables; Standard errors are clustered at the LMR-level. Source: SOEP v.37, 2009-2020."),
       file = paste0(tables_path, "robustness_heterogeneity/Appendix_Table_D29b.tex"))



#-------------------------------------------------------------------------------# 
#-------------------------------------------------------------------------------# 
# Clear workspace
rm(list=ls())
#-------------------------------------------------------------------------------# 
#-------------------------------------------------------------------------------# 



